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Foreword 

The  goal  of  this  project  was  to  address  basic  nonlinear  problems  in  the  modeling  of  physical 
systems  with  engineering  relevance.  Nonlinear  conservation  laws  were  chosen  for  study  because 
these  differential  equations  are  basic  to  many  applications.  We  also  studied  stochastic  issues,  as 
they  are  often  of  great  importance  in  nonlinear  modeling  of  physical  phenomena.  In  these  studies, 
we  used  mathematical  theory,  first  principle  computations  and  analysis  of  experimental  and  field 
data. 

Motivating  applications  of  our  work  included,  inertial  confinement  fusion,  shock  and  gravity 
induced  fluid  mixing,  fluid  jet  breakup,  ground  water  flows,  petroleum  reservoir  modeling,  and 
high  velocity  impacts.  Additional  applications  included  vehicle  signature  identification  (radar  cross 
sections),  advanced  manufacturing  technology,  pharmacology  and  drug  activity  through  protein- 
DNA  binding  studies,  and  remediation  of  environmentally  polluted  facilities. 
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3.  List  of  Appendixes,  Illustrations,  and  Tables 

The  tracked  fronts  for  the  height  of  burst  problem  with  a  thermal  boundary  layer. 

(a)  shows  the  entire  grid,  and  (b)  a  blow  up  of  the  lower  right  portion  of  the  grid. 

The  vertical  coordinates  have  been  stretched  by  a  factor  of  67%  in  (b).  The  grid  is 
50  X  50,  so  that  the  boundary  layer  is  exactly  five  mesh  zones  in  width.  The  frame 

(b)  occupies  a  6  x  6  grid  subregion  of  the  full  computation .  5 

Graphs  of  the  ramp  angle  vs.  the  trajectory  angle  of  the  Mach  node,  0^  is  the 
detachment  angle,  beyond  which  regular  reflection  is  theoretically  impossible,  and 
the  mechanical  equilibrium  or  von  Neumann  point.  The  experimental  results  can 
get  no  closer  to  transition  than  seven  degrees,  and  the  CRAY  results  no  closer  than 
three  degrees,  both  due  to  a  lack  of  resolution  near  the  wall.  The  front  tracking 
results  represent  resolution  to  within  one  half  degree  of  transition .  6 

(a)  The  amplitude  growth  rate,  a(^),  of  a  shocked  air-SFe  interface.  This  graph 
compares  the  results  of  experimental  averages,  front  tracking  simulation,  linear 
theory  and  Richtmyer’s  impulsive  model.  Also  shown  are  results  of  a  least  squares 
fit  to  the  front  tracking  amplitude  data  over  the  period  of  experimental  observation, 
(b-e)  Pressure  plots  at  a  series  of  times  for  the  nonlinear  air-SFa  solution.  A  cascade 
of  shock  waves  generated  by  the  self-interaction  of  the  transmitted  and  reflected 
waves  propagate  back  toward  the  interface  and  affect  the  perturbation  growth  rate 
at  early  and  intermediate  times .  8 

In  the  experiments  the  mold  was  filled  till  “breakthrough”,  after  which  time  the 
resin  cured  and  the  void  distribution  was  measured  along  the  mold.  ..The  voids 
concentrate  on  the  outlet  side  of  the  mold.  Our  relative  permeability  models  have  a 
few  physical  parameters.  The  figure  shows  that  these  parameters  can  be  adjusted  to 


give  quantitative  agreement  between  experiment  and  simulation .  10 

Comparison  of  the  experimental  data  and  our  calculations  with  the  new  material 
model  with  the  new  material  parameter  C2  =  0.035  Mb  ms.  The  plots  display  the 
velocity  of  the  right-hand  free  surface  vs.  time  in  tests  161,  232,  and  390 .  12 
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4.  Scientific  Activities 

L  Research  Summary 

LI.  Global  Solutions  to  the  Compressible  Euler  Equations  with  Geometrical  Struc¬ 
ture.  Joint  work  with  G.-Q.  Chen  [Chen  and  Glimm,  1994a,  Chen  and  Glimm,  1994b]  settled  a 
theoretical  question  that  has  been  outstanding  for  several  decades,  namely  how  to  control  the 
influence  of  reflected  waves  from  infinity,  to  prove  existence  theorems  for  the  Euler  equations  of 
isentropic  gas  dynamics  for  radially  symmetric  flow  in  two  and  three  dimensions.  The  domain 
considered,  r  >  1,  includes  the  singularities  associated  with  infinity,  but  excludes  the  singularity 
of  focusing  at  the  origin.  The  main  idea  in  the  proof,  which  uses  the  compensated  compactness 
framework,  is  to  construct  approximate  solution  operators  and  approximate  Riemann  solutions 
that  incorporate  leading  order  radial  corrections.  These  are  based  on  exact  solutions  of  the  radial 
equations  which  are  steady,  i.e.  no  with  no  time  dependence. 

1,2.  Stochastic  Analysis  of  Fluid  Waves.  Macroscopic  behavior  is  the  result  of  averages  over 
small  scale  fluctuations.  Ensemble  averages  are  a  convenient  formalism  to  express  this  principle. 
Ensemble  averages  introduce  the  severe  complication  of  function  space  integrals,  or  of  analysis  in 
infinite  dimensional  function  spaces.  Such  difficulties  are  also  encountered  in  statistical  physics 
and  quantum  field  theory  [Glimm  and  Jaife,  1981],  We  used  four  methods  to  study  ensemble  av¬ 
erages:  renormalization  group  methods,  moment  expansions,  perturbation  theory  (ordinary  and 
renormalized),  and  direct  numerical  simulation. 

1.2.1.  Renormalization  Group  Analysis.  In  the  renormalization  group  (RNG)  approach,  there  is  a 
continuum  of  length  scales,  with  integration  over  all  small  scales  used  to  define  the  effective  equa¬ 
tions  on  some  intermediate  length  scale.  The  introduction  of  new  (small)  length  scales  and  the 
change  in  location  of  the  intermediate  length  scale  is  treated  differentially,  to  give  the  renormahza- 
tion  group  equation.  The  case  of  a  fixed  point  for  this  equation  is  of  special  interest.  In  this  case, 
the  problem  formulation  and  equations  do  not  change  as  the  intermediate  length  scale  is  varied. 
The  equations  and  solutions  are  self  similar.  Power  law  exponents  and  premultiplying  coefficients 
then  determine  the  solution.  Usually  the  exact  integration  of  the  small  length  scales  gives  rise  to 
extremely  complicated  effective  equations,  which  are  difficult  to  work  with.  Fortunately,  in  the  case 
of  a  fixed  point,  approximations  will  not  change  the  exact  answer.  The  fixed  point  will  normally 
be  stable  to  perturbations,  or  at  least  to  many  of  them.  Let  us  call  the  stable  perturbing  directions 
for  the  equations  of  motion  irrelevant,  and  the  unstable  directions,  relevant.  Then  any  irrelevant 
terms  in  the  exactly  integrated  effective  equations  can  be  discarded.  The  effective  equations  at  the 
intermediate  length  scale  will  contain  only  relevant  terms.  The  three  step  renormalization  group 
process  is  as  follows:  integrate  a  (differential)  unit  of  length  scale  dk,  rescale,  so  that  the  equations 
are  formulated  at  a  constant  length  scale,  and  finally  truncate,  so  that  irrelevant  terms  are  dropped. 
If  this  process  yields  the  identity,  so  that  the  equations  do  not  change  under  integration  of  length 
scales,  then  a  fixed  point  ha.s  been  found. 

We  have  use  a  variant  of  this  method  to  study  the  Rayleigh- Taylor  fluid  mixing  process.  This 
mixing  process  is  an  acceleration  driven  instability  of  an  interface  separating  fluids  of  different 
densities.  We  wish  to  characterize  the  height 

h  =  aAgt^ 

of  the  mixing  zone  as  a  function  of  time.  L  where  g  is  the  (gravitational)  acceleration,  A  ~  {p2  — 
Pi)/{p2  +  Pi)  is  the  Atwood  number,  a  buoyancy  renormalization  of  g,  the  pi  are  densities  on  the 
two  sides  of  the  interface,  and  a  is  a  dimensionless  constant  that  characterizes  the  mixing  rate.  For 
incompressible  Hows,  a  =  .06  is  a  universal  constant.  An  unusual  feature  of  this  problem  is  that 
the  time  evolution  enforces  dynamically  a  change  (increase)  of  the  effective  length  scale.  Thus  time 
can  be  taken  as  the  otherwise  artificial  parameter  in  the  renormalization  group  equation  dynamics. 
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The  temporal  dynamics  coincides  with  the  renormalization  group  dynamics  and  the  renormalization 
group  describes  the  transients,  z.e.  the  dynamical  approach  to  steady  state  self  similarity.  We  feel 
that  this  is  an  important  observation,  that  may  find  validity  in  many  related  “totally  unstable, 
unstable  on  all  length  scales”  instabilities,  such  as  the  Taylor-SaflFman,  KelvimHelmholtz,  and 
metastable  phase  transition  instabilities. 

Our  RNG  study  of  the  Rayleigh-Taylor  instability  starts  with  a  modification  of  the  Sharp- 
Wheeler  model,  and  defines  a  statistical  model  for  the  dynamic  evolution  of  an  ensemble  of  bubbles. 
This  model  is  regarded  as  an  approximation  to  the  two  fluid  Euler  equations,  describing  the  fluid 
evolution  directly.  The  bubble  model  has  its  own  laws  of  dynamical  evolution.  Each  bubble  moves 
as  the  sum  of  two  velocities.  The  first  is  due  to  the  bubble  itself,  and  is  the  terminal  velocity  taken 
from  single  bubble  studies  of  the  Rayleigh-Taylor  problem.  It  is  a  function  of  the  bubble  width. 
The  second  contribution  is  a  velocity  associated  with  the  bubble  interactions,  and  is  defined  in 
terms  of  the  relative  position  of  the  bubble  and  its  neighbor.  It  is  also  derived  from  single  bubble 
studies  of  the  Rayleigh-Taylor  problem.  At  sufficiently  large  relative  heights,  a  merger  process 
will  occur  between  adjacent  bubbles,  leading  to  a  new  bubble,  of  width  equal  to  the  combined 
widths  of  the  two  that  are  merging.  The  parameter  in  the  model  that  controls  merger  is  set  from 
numerical  Rayleigh-Taylor  studies  of  two  adjacent  bubbles,  and  checked  by  analysis  of  experimental 
data.  This  model  [Glimm  and  Sharp,  1990,  Glimm  et  ah',  1991],  is  then  approximated.  First,  the 
adjacent  bubble  interactions  are  approximated  by  random  pairs  selected  from  the  ensemble,  and 
then  the  probability  distributions  for  the  bubbles  are  approximated  by  a  specific  parametric  form. 
The  dynamics  is  then  projected  into  a  dynamics  for  the  parameters  of  the  bubble  probabilities. 
This  dynamics  is  the  renormalization  group  equation,  and  a  fixed  point  is  found.  Further  analysis 
showed  that  this  fixed  point  is  in  quantitative  agreement  with  the  growth  rate  of  the  mixing  zone,  as 
determined  experimentally  [Read,  1984,  Youngs,  1984a],  and  computationally,  by  the  front  tracking 
method  [Zhang,  1990]. 

In  this  application  of  the  renormalization  group  method  to  the  Rayleigh-Taylor  problem,  we  did 
not  analyze  relevant  and  irrelevant  terms.  Rather  we  followed  the  methods  of  intuitive  physics 
to  justify  the  approximations.  Our  approach  has  the  benefit  that  each  intermediate  step  has  a 
physical  interpretation.  To  compensate  for  the  omission,  an  independent  check  on  the  conclusions 
was  required,  which  we  supplied  in  the  form  of  agreement  of  the  RNG  fixed  point  behavior  with 
experimental  and  numerical  simulation  data. 

1.2.2.  Moment  Expansions.  In  moment  expansion  methods,  dynamical  equations  are  derived  for 
low  order  moments  of  the  ensemble  average  [Drew,  1983,  McComb,  1990].  Equation  nonlinearities 
couple  low  order  moments  to  higher  ones,  and  a  closure  hypothesis  is  needed  to  complete  the 
dynamical  equations.  This  hypothesis  is  not  to  be  “derived”  in  a  mathematical  manner  from 
existing  laws  of  physics,  but  is  to  be  “introduced”  as  a  physical  principle  on  its  own  merit.  Thus 
it  must  be  justified  by  appeal  to  measured  or  simulated  data,  and  its  use  must  be  restricted  to 
regimes  for  which  it  is  valid.  Usually  the  limitations  of  the  moment  expansion  method  lie  in  this 
step  of  validation  and  restriction. 

We  have  applied  the  moment  expansion  method  to  the  study  of  the  Rayleigh-Taylor  mixing  layer. 
Our  preliminary  results  are: 

(a)  A  two  phase  description  fits  the  data  better  than  a  turbulence  description. 

(b)  A  new  first  moment  two  phase  closure  is  proposed. 

(c)  A  one  parameter  family  of  solutions  is  obtained  to  the  ensemble  averaged  equations. 

(d)  A  better  understanding  of  the  loss  of  universality  for  the  compressible  Rayleigh-Taylor 
mixing  phenomena  is  achieved. 

(e)  The  new  equations  are  hyperbolic,  z.e.,  they  have  purely  real  characteristics. 

We  also  note  earlier  work  [Stewart  and  WendroAF,  1984]  in  which  purely  hyperbolic  multiphase 
equations  are  proposed. 
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The  new  two  phase  ensemble  averaged  equations  have  the  form: 

dai,  _  _  ,dak  ^ 

—  i-  {aiV2  4-  a2Vi)-^  =  0  , 


dO'kPk  ^^kPk'^k 
dt  dz 

dakpkVk  diakpkVkVk) 


-0  , 
d{akpk) 


dt 


dz 

-(arp2  +  a2P\) 


dz 


+  <^kPk9 


d{akpkh)  ,  d(o^kPkhvk) 

+  - -  ==  ~Pk 


dak 
dz 

-  d{akVk) 


for  A:  =  1,  2  and 


dt 


dz 


dz 

O:'!  +  0^2  —  1  5 


+  {aiP2V2  +  a2PiVi)^§^ 
dz 


where  is  the  volume  fraction,  bars  denote  volume  weighted  averages,  and  tildes  denote  density 
weighted  averages,  p  is  density,  v  is  velocity,  p  is  pressure,  g  gravity  and  e  is  the  internal  energy.  See 
[Harlow  and  Amsden,  1975,  Youngs,  1984b.  Freed  et  al,  1991]  for  a  discussion  of  equations  (deriv¬ 
able  from  the  above)  in  which  there  is  a  single  equilibrated  pressure,  and  the  flow  is  incompressible. 

We  studied  both  turbulence  and  two  phase  descriptions  of  this  mixing  data,  as  both  have  been 
proposed  for  this  problem.  We  derived  an  identity  to  write  the  second  order  turbulence  moments 
such  as  the  Reynolds  stress  as  products  of  the  first  order  two  phase  moments,  plus  an  intrinsic  two 
phase  second  order  turbulence  contribution.  The  latter  was  seen  to  be  negligible  in  the  Rayleigh- 
Taylor  data  we  analyzed.  This  fact  supports  our  conclusion  that  the  two  phase  description  is  more 
appropriate  for  this  data  [Chen  et  ah,  1995]. 


1.2.3.  Perturbation  Theory.  The  method  of  perturbation  theory  is  useful  for  small. random  fluctu¬ 
ations.  Together  with  coworkers,  we  applied  both  ordinary  and  renormalized  perturbation  theory 
to  the  problem  of  transport  by  a  random  velocity  field. 

The  transport  of  a  passive  (noninteracting)  scalar  contaminant  u  by  a  random  velocity  field  v 
arises  in  environmental  studies  of  ground  water.  The  equation  governing  this  flow  [Dagan,  1989]  is 

Uf  A-  V  ■  vu  ~  0  . 


Here  if  is  a  Gaussian  random  field,  defined  in  terms  of  its  covariance,  {v{x).,v{y)).  We  assume  that 
the  statistics  are  stationary,  self  similar  (fractal),  isotropic  and  incompressible.  More  generally, 
random  fields  with  stationary  increments  [Gelfand  and  Vilenkin,  1964,  Yaglom,  1986]  have  been 
used  to  model  geostatistics;  in  some  approximation,  these  random  flelds  describe  the  velocities  as 
well.  It  follows  that  the  covariance  has  a  povrer  law  form,  as  a  function  of  r  =  1^;  —  yj,  namely  br~^. 
Since  u  is  a  random  field,  so  is  li  =  u{x,t).  and  the  problem  is  to  relate  the  statistics  of  v  to  that 
of  u.  The  stochastic  property  of  v  reflects  the  lack  of  knowledge  of  details  of  small  scale  variations 
in  the  geology,  and  hence  in  the  ground  water  flow  velocity,  v.  We  write  u  =  u  -f  where  v  =  {v) 
is  the  ensemble  average,  or  mean  flow.  Perturbation  theory  expresses  u  as  a  series  in  powers  of 
6v.  The  leading  order  contribution  to  w  resulting  from  statistical  fluctuations  Sv  in  v  is  a  diffusion 
term,  added  to  the  transport  of  u, 

Ut  V  '  vu  ~  u/\u  , 


where 


u  ~ 


J {6v{x  +  vt)6v{x))dt  , 


In  the  case  of  slowly  decaying  velocity  correlations  the  dispersion  is  anomalous,  i.e.  not  Fickean, 
and  u  is  a  solution  of  the  diffusion  equation  with  a  time  dependent  diffusion  constant.  Anomalous 
diffusion  is  consistent  with  field  data,  for  ground  water  flow,  for  which  it  is  known  that  the  disper¬ 
sion  coefficients  increase  with  time  or  travel  distance  [Gelhar  et  al.,  1985],  by  up  to  six  orders  of 
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magnitude.  Our  analysis  traced  this  anomalous  aspect  of  the  dispersion  in  ground  water  flows  to 
slowly  decaying  correlations  in  the  geological  data,  specifically  to  the  correlations  of  permeability, 
and  hence  to  the  velocity  field,  in  our  models.  Our  main  conclusion  is  that  a  self  similar,  or  fractal 
description  of  the  velocity  field  leads  to  mixing  behavior  that  is  not  self  similar,  but  whose  short 
and  long  distance  limits  have  (different,  but  explicitly  determined)  scaling  laws.  Let  the  mixing 
length  I  be  a  function  of  the  travel  distance  L.  An  asymptotic  scaling  law  states  that  I  ~  L^.  The 
scaling  law  a  =  1/2  is  Fickean,  and  is  valid  for  large  L  if  /3  >  1,  i.e.  if  the  integral  defining  v  above 
is  convergent.  Otherwise,  the  diffusion  is  not  Fickean.  and,  in  general, 


For  ground  water  problems,  the  transient  i.e.  preasymptotic  dispersion  is  more  important  than 
the  asymptotic  dispersion.  See  [Furtado  et  ah,  1992,  Glimm  et  ah,  1992,  Glimm  et  ah,  1993a]. 

1.3.  Computation.  Recent  developments  have  enabled  our  front  tracking  code  to  compute  suc¬ 
cessfully  several  problems  of  long  standing,  leading  to  greatly  improved  computations  and,  for  some 
cases,  a  first  correct  computation.  These  new  capabilities  are  also  laying  the  groundwork  for  new 
fundamental  studies  in  the  future. 

1.3.1.  Stochastic  Flows.  The  method  of  direct  numerical  simulation  is  a  powerful  tool  for  the  anal¬ 
ysis  of  stochastic  theories.  It  provides  data  for  assessing  the  validity  of  approximate  theoretical 
models  and  it  provides  solutions  for  which  there  is  no  applicable  theory.  We  have  studied  the 
two  problems  introduced  above  by  this  method:  the  Rayleigh- Taylor  problem  and  the  anomalous 
diffusion  due  to  transport  by  a  random  velocity  field.  Both  studies  focused  on  the  growth  rate  of 
the  mixing  zone.  For  the  Rayleigh-Taylor  problem,  we  have  the  only  compressible  computation 
which  shows  agreement  with  incompressible  experiments  in  the  incompressible  limit.  (There  is  no 
compressible  experimental  data  available.)  We  discovered,  on  the  basis  of  computational  experi¬ 
ments,  a  loss  of  universality  for  the  growth  rate  in  the  compressible  case,  and  a  possible  doubling 
of  the  growth  rate  [Chen  et  ah,  1993,  Deng  et  ah,  1993].  In  more  recent  work,  we  obtained  an  un¬ 
derstanding  of  the  loss  of  universality,  or  ensemble  dependence  of  the  growth  rate.  It  appears  that 
the  initial  instability  amplitude,  or  the  relative  amplitudes  of  the  long  and  short  wave  lengths  in 
the  initial  data  are  significant  variables  for  the  growth  rate  determination  [Chen  et  ah,  1995]. 

Our  simulations  displayed  anomalous  dispersion  in  agreement  with  both  ordinary  and  renormal¬ 
ized  perturbation  theory  solutions,  for  transport  by  a  random  velocity  field  [Furtado  et  ah,  1992, 
Glimm  et  ah,  1993a]. 

In  [Li  et  ah.  1995]  we  compared  front  tracking  to  TVD  schemes,  with  and  without  artificial 
compression.  The  latter  was  also  discussed  in  three  dimensions.  It  was  found  that  TVD  is  diffusive 
for  the  Rayleigh-Taylor  problems,  and  that  this  problem  was  to  a  large  extent  overcome  by  the  use 
of  artificial  compression.  The  level  set  method,  also  used  in  this  study,  is  a  post  processing  step, 
and  affects  the  graphical  output,  but  not  the  computation. 

1.3.2.  Shock  Diffraction.  The  first  of  these  new  developments  is  the  use  of  front  tracking  for  shock 
wave  diffraction  by  a  fluid  interface  or  by  an  obstacle.  The  novelty  is  not  the  principle  of  such  a 
front  tracking  computation,  which  was  first  achieved  by  the  proposers  and  coworkers  some  years 
ago,  but  its  robustness.  Robustness  is  measured  by  the  number  of  cases  handled  automatically, 
or  by  the  complexity  of  the  problems  thus  solved  automatically.  The  computations  are  based  on 
shock  polars.  which  are  equations  defined  by  shock  wave  jump  conditions,  composed  in  sequence 
for  all  the  waves  meeting  at  a  point.  As  a  tool  for  graphical  understanding,  the  shock  polars  are 
projected  on  some  convenient  plane  (such  as  the  pressure,  turning  angle  plane).  We  use  shock  polar 
analysis  as  a  constructive  numerical  step  in  our  front  tracking  finite  difference  code.  The  analysis 
inserts  analytic  information  into  the  computation  where  it  is  most  needed,  at  the  points  of  greatest 
computational  difficulty.  The  result  of  this  method  are  computations  of  remarkable  accuracy,  often 


ADV,  COMP.  Sz  MODELING  FOR  NONLINEAR  STOCHASTIC  WAVES 


5 


on  very  coarse  ^^rids.  The  dynamic  change  of  the  topology  defined  by  the  fronts  is  automated  to 
allow  the  formation  of  new  shock  waves,  and  the  elimination  of  others  as  they  disappear  from  the 
computation  or  become  less  important.  The  hierarchy  of  possible  shock  diffraction  diagrams  is 
exceedingly  complex,  and  not  all  of  them  are  presently  supported  as  front  tracking  capabilities. 
Rather,  support  for  enough  of  them  has  been  achieved  to  allow  solution  of  various  interesting 
problems.  Validation  studies  for  a  regular  Mach  reflection,  computed  by  these  methods,  have 
been  reported  [Boston  et  ah,  1993b].  One  interesting  flow  pattern  studied  was  the  height  of  burst 
problem  for  the  interaction  of  an  expanding  spherical  shock  wave  with  a  thermal  boundary  layer 
above  a  planar  rigid  surface,  see  Figure  1.1. 
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Figure  1.1.  The  tracked  fronts  for  the  height  of  burst  problem  with  a  thermal 
boundary  layer,  (a)  shows  the  entire  grid,  and  (b)  a  blow  up  of  the  lower  right 
portion  of  the  grid.  The  vertical  coordinates  have  been  stretched  by  a  factor  of  67% 
in  (b).  The  grid  is  50  x  50,  so  that  the  boundary  layer  is  exactly  five  mesh  zones 
in  width.  The  frame  (b)  occupies  a  6  x  6  grid  subregion  of  the  full  computation. 

Another  interesting  study  concerned  the  regular  reflection  -  Mach  stem  transition.  Computations 
with  wedge  angles  within  a  fraction  of  a  degree  of  the  von  Neumann  transition  angle  were  obtained. 
The  best  comparison  computation  (based  on  Godunov  methods,  local  mesh  refinement  and  a  limited 
version  of  front  tracking)  [Henderson  et  ah.  1993]  required  considerably  more  resources  to  obtain 
a  less  precise  determination  of  the  transition  point.  The  two  computations  agreed  within  their 
common  domain  of  definition.  See  Figure  1.2. 

The  most  important  application  of  this  enhanced  front  tracking  shock  wave  diffraction  capability 
is  to  shock  induced  turbulent  mixing,  described  below. 

1.3.3.  Parallelism.  The  second  recent  development  for  front  tracking  is  its  extension  to  parallel 
MIMD  hardware.  This  extension  has  enabled  most  of  the  other  computations  reported  here.  For 
hyperbolic  conservation  laws,  with  explicit  time  step  integration,  domain  decomposition  is  imple¬ 
mented  with  use  of  ‘‘ghost  cells”  to  extend  the  size  of  each  domain  by  approximately  one-half  of  the 
stencil  width.  The  data  from  the  ghost  cells  allow  the  ordinary  serial  code  to  compute  for  a  time 
step  on  each  domain.  At  the  end  of  the  time  step,  communication  is  used  to  update  the  ghost  cells, 
as  a  block.  In  this  method,  communication  is  restricted  to  data  proportional  to  the  boundary  of  the 
domains,  and  the  number  of  communication  packets  is  0(1).  A  complication  of  front  tracking,  for 
this  algorithm,  is  its  extensive  use  of  the  C  language  data  structures,  including  function  and  data 
pointers.  After  communication,  the  data  pointers  are  incorrectly  set,  and  must  be  reset  through 
controlled  storage  allocation  within  designated  storage  blocks  and  through  pointer  arithmetic,  to 
generate  the  oflset  in  relative  addresses  between  distinct  blocks.  In  addition,  a  number  of  previ¬ 
ously  solved  problems,  such  as  the  remeshing  of  points  on  the  front,  had  to  be  reexamined,  since 
the  previously  used  global  algorithm  of  uniform  spacing  by  arclength  is  not  convenient  relative 
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Figure  1.2.  Graphs  of  the  ramp  angle  9-^^^  vs.  the  trajectory  angle  of  the  Mach 
node,  X-  is  the  detachment  angle,  beyond  which  regular  reflection  is  theoretically 
impossible,  and  6f^  the  mechanical  equilibrium  or  von  Neumann  point.  The  exper¬ 
imental  results  can  get  no  closer  to  transition  than  seven  degrees,  and  the  CRAY 
results  no  closer  than  three  degrees,  both  due  to  a  lack  of  resolution  near  the  wall. 

The  front  tracking  results  represent  resolution  to  within  one  half  degree  of  transition. 

to  one  which  is  local  to  each  processor.  The  parallelization  has  been  a  success  by  all  measures: 
The  parallel  speed  up  is  above  90%.  The  time  spent  in  communications  is  small  relative  to  the 
time  spent  in  computations.  The  results  produced  are  of  considerable  scientific  interest,  and  could 
not  have  been  obtained  without  the  parallelization.  The  level  of  utilization  of  our  local  parallel 
hardware  approaches  100%. 

1.3,4.  Three  Dimensional  Front  Tracking,  The  third  recent  development  is  the  extension  of  front 
tracking  to  three  dimensions.  This  extension  is  still  undergoing  validation  studies,  and  a  paper 
describing  its  capabilities  is  being  prepared.  See  also  [Li  et  ah,  1993,  Glimm  et  ah,  1995a].  We 
discuss  in  detail  here  one  subroutine  of  this  program,  namely  the  efficient  interpolation  of  piecewise 
smooth  data  in  three  dimensions.  The  front  tracking  interface  structure  will  support  the  description 
of  an  arbitrary  set  of  discontinuity  surfaces,  intersecting  along  curves,  which  themselves  intersect 
on  codimension  three  objects,  i.e.  points,  and  referred  to  here  as  nodes.  The  interface  structure 
is  represented  discretely  as  piecewise  linear  objects.  Thus  the  discretized  surface  is  composed  of 
triangles,  and  the  curves  are  composed  of  line  segments,  called  bonds.  The  two  ends  of  the  line 
segments  (bonds)  are  called  points.  The  vertices  of  the  surface  triangles  are  also  called  points,  and 
for  surfaces  meeting  in  a  common  curve,  there  will  also  be  triangles  meeting  at  a  common  bond. 
Geometrically  identical  points  are  shared  as  a  common  object.  The  surface  is  oriented,  so  that 
a  distinguished  normal  direction  (positive)  is  introduced.  At  each  point  there  is  a  positive  and  a 
negative  storage  location  for  discontinuous  solution  state  variables.  Along  curves,  where  multiple 
surfaces  meet,  states  are  stored  for  positive  and  negative  sides  of  each  surface  bounded  by  the  curve, 
at  each  point  (bond  start  or  endpoint)  of  the  curve.  In  addition,  state  variable  storage  is  provided 
at  all  regular  grid  locations,  at  grid  cell  centers.  The  interpolation  problem  is  to  construct  an 
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efficient  solution  function  that  will  return  an  interpolated  solution  at  an  arbitrary  location  in  three 
dimensional  space,  and  that  is  guaranteed  robust,  in  the  following  two  senses:  The  discontinuity 
interface  will  divide  space  into  connected  components.  If  the  connected  component  is  specified, 
then  the  interpolation  is  guaranteed  to  be  among  grid  cell  values  associated  with  this  component 
only.  If  no  component  is  given,  then  the  component  associated  with  this  location  will  be  computed 
efficiently,  and  the  previous  method  used. 

Efficient  precomputation  of  components  and  their  storage  on  those  grid  cells  not  meeting  the 
interface  (so  that  they  are  uniquely  defined)  is  accomplished  by  a  hashing  method,  which  lists  all 
triangles  meeting  a  given  mesh  cell.  From  this  list,  an  0(1)  computation  will  give  local  component 
information.  This  local  information  can  be  propagated  to  nearest  neighbors,  if  there  is  no  interven¬ 
ing  interface,  and  thus  extended  to  all  of  the  computational  domain.  The  evaluation  of  the  robust 
solution  function  uses  the  dual  grid,  with  data  located  at  grid  corners,  not  at  grid  cell  centers.  Eval¬ 
uation  is  by  bilinear  interpolation  for  points  lying  in  regular  dual  grid  cells.  For  grid  cells  meeting 
the  interface,  a  further  tetrahedralization  is  needed.  Because  the  tetrahedra  will  be  used  for  linear 
interpolation,  no  new  interior  vertices  (at  which  there  would  be  missing  data)  can  be  introduced.  It 
is  known  that  the  above  problem  is  in  general  insoluble,  unless  new  points  are  introduced.  In  fact, 
it  would  be  possible  to  introduce  new  points  on  the  interface,  giving  it  very  small  triangles.  Data 
at  the  new  triangular  points  would  be  obtained  by  linear  interpolation  on  the  original  triangles. 
This  solution  is  not  very  desirable,  as  it  appears  to  be  computationally  ineflScient.  Instead,  we  use 
a  Dulaney  tetrahedralization  of  space.  This  tetrahedralization  will  not  respect  the  division  of  space 
into  components  by  the  interface.  In  other  words,  an  individual  tetrahedron,  with  vertices  defined 
by  one  component,  may  intersect  another  component.  Thus  we  divide  our  data  points  into  the 
components  they  represent,  and  construct  an  independent  tetrahedralization  for  each  component. 
For  efficiency  reasons  this  is  done  locally  within  each  dual  grid  mesh  block.  Since  the  interpolating 
function  is  to  be  supplied  with  a  component  as  an  argument,  or  by  a  side  computation,  it  will  call 
up  the  interpolation  based  on  the  Dulaney  tetrahedralization  of  that  component. 

1.3.5.  E&M  Scattering  Cross  Sections.  We  also  developed  a  two  step  multigrid  algorithm  for  the 
solution  of  electro-magnetic  scattering  problems  by  boundary  integral  methods.  These  problems 
are  memory  limited,  and  since  only  the  coarse  grid  inverse  matrix  is  stored  (the  fine  grid  direct 
matrix  is  computed  on  the  fly,  using  tabulated,  efficient  special  functions),  problems  up  to  25  times 
larger  than  normal  (single  step)  methods  can  be  solved  by  this  algorithm.  We  also  investigated 
limitations  of  this  algorithm,  associated  with  resonances.  For  dielectric  material,  and  also  for 
indentations  in  a  perfectly  conducting  plane  surface,  the  problem  of  resonances  was  overcome.  The 
algorithm  was  implemented  as  MIMD  parallelized  code.  On  a  small  parallel  computer,  we  are 
able  to  solve  30.000  equations  using  a  two  level  iterative  approach,  for  dielectric  materials  and  for 
indentations  in  a  plane  surface  [Asvestas  et  al.,  1993,  Asvestas  et  ah,  1994].  This  same  code  has 
been  adapted  to  solve  radiative  heat  transport  problems,  for  application  to  crystal  growth  modeling. 
Heat  transport  is  of  interest  to  us  as  a  detail  in  a  planned  integrated  approach  to  the  modeling, 
design,  and  manufacturing  control  of  the  crystal  pulling  process, 

1.4.  Applications. 

1.4.1,  Shock  Induced  Mixing.  The  proposers  and  coworkers  have  recently  obtained  significant  progress 
in  understanding  turbulent  mixing  and  multiphase  flow.  Our  efforts  have  been  concentrated  on  ac¬ 
celeration  induced  instabilities  of  an  interface  between  fluids  of  differing  densities.  The  case  of 
impulsive  acceleration,  as  due  to  a  shock  wave,  is  known  as  the  Richtmyer- Meshkov  instability; 
the  case  of  steady  acceleration,  as  due  to  a  gravitational  force,  is  known  as  the  Rayleigh-Taylor 
instability. 

The  agreement  obtained  (for  the  first  time)  between  computation  and  experiment  in  the  growth 
rate  of  the  Richtmyer- Meshkov  instability  [Grove  et  al.,  1993a,  Grove  et  al.,  1993c,  Grove,  1994, 
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Grove  et  al.,  1993c.  Boston  et  al.,  19930,  Boston  et  al.,  1993a]  is  the  most  important  aspect  of  the 
progress  reported  here.  In  Figure  L3a  we  show  the  growth  rates  for  a  Richtmyer-Meshkov  insta¬ 
bility,  determined  by  a  front  tracking  simulation.  Superimposed  on  this  plot  is  the  growth  rate 
reported  from  experiment  [Benjamin  et  ah,  1993],  that  obtained  by  linear  perturbation  theory,  and 
the  growth  rate  of  the  Richtmyer  impulsive  model. 

These  computations  have  produced  a  far  deeper  understanding  of  the  theory  of  this  instability 
than  has  been  previously  available  [Boston  et  al..  1995a,  Grove  et  ah,  1995,  Holmes  et  ah,  1995, 
Boston  et  ah.  1995b].  We  understand  the  shock  diffraction  mechanisms  that  yield  the  decay  in  the 
mixing  growth  rates.  In  Figure  1.3b-e,  we  show  plots  of  pressure  at  the  times  corresponding  to  the 
down  turns  in  the  growth  rate,  as  shown  in  Figure  1.3a.  Note  the  appearance  of  a  pressure  pulse 
at  either  the  peak  or  the  trough  of  the  mixing  interface  in  these  figures,  to  cause  a  sudden  dip  in 
the  growth  of  the  peak  to  trough  separation.  The  pressure  pulses  originate  in  the  intersection  of 
curved  shock  waves,  moving  transversely  to  the  interface  and  coming  from  the  curvature  of  the 
original  interface. 


Figure  1.3.  (a)  The  amplitude  growth  rate,  a(i),  of  a  shocked  air-SFe  interface. 
This  graph  compares  the  results  of  experimental  averages,  front  tracking  simulation, 
linear  theory  and  Richtmyer's  impulsive  model.  Also  shown  are  results  of  a  least 
squares  fit  to  the  front  tracking  amplitude  data  over  the  period  of  experimental 
observation,  (b-e)  Pressure  plots  at  a  series  of  times  for  the  nonlinear  air-SFe  solu¬ 
tion.  A  cascade  of  shock  waves  generated  by  the  self-interaction  of  the  transmitted 
and  reflected  waves  propagate  back  toward  the  interface  and  affect  the  perturbation 
growth  rate  at  early  and  intermediate  times. 


We  also  understand  the  failure  of  other  computations  to  give  the  right  answer:  numerical  diffusion 
of  the  interface  led  to  a  diminished  shock-interface  interaction.  The  numerical  diffusion  becomes 
more  pronounced  as  the  computation  progresses,  so  that  the  effects  of  moderate  or  late  time  shock- 
interface  interactions  were  especially  diminished.  However,  it  was  these  late  time  interactions,  due 
to  reflections  from  the  curved  contact  discontinuity  and  the  curved  shock  waves,  that  caused  the 
decay  in  the  growth  rate  of  the  instability.  Lagrangian  computations  [Meyer  and  Blewett,  1972], 
with  a  sharp  fluid  interface,  were  not  carried  to  late  time,  presumably  due  to  mesh  tangling  prob¬ 
lems.  An  incorrect  physical  picture,  corrected  in  work  of  the  proposers  and  coworkers  cited  above, 
held  that  the  growth  rate  for  the  instability  was  constant  in  time.  This  picture  led  to  inappropriate 
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comparisons  between  growth  rates,  for  example  between  those  observed  at  different  time  intervals, 
and  it  led  to  the  recording  of  a  single  growth  rate  in  experiment,  rather  than  the  recording  of  the 
growth  history  over  the  longest  possible  observable  interval. 

Other  computations  had  reasonable  performance  for  early  time  growth  rates,  but  they  developed 
inaccuracies  before  the  period  of  experimental  observations  began  and  displayed  errors  in  over 
predicting  the  growth  rate  during  the  experimentally  observed  period,  typically  by  a  factor  of  two. 
Popular  theoretical  models  for  prediction  of  this  instability  also  give  incorrect  answers,  and  the 
basis  for  these  discrepancies  is  now  understood.  Fundamentally,  these  methods  are  based  on  small 
amplitude  assumptions  and  linearization  of  the  perturbation  of  an  idealized  flat  interface.  This 
assumption  is  invalid  for  the  times  and  instability  amplitudes  considered, 

1.4.2.  Resin  Transfer  Molding.  Resin  Transfer  Molding  (RTM)  is  a  process  to  produce  light  weight, 
high  strength  structural  components,  possibly  with  complex  geometries,  for  aerospace  and  auto¬ 
motive  applications.  In  this  process,  a  liquid  plastic  resin  is  injected  under  pressure  into  a  woven 
mat  of  carbon  fiber;  the  part  is  then  set.  or  hardened,  by  curing,  i.e,  polymerization  of  the  resin, 
by  baking  it  in  an  autoclave.  There  are  many  variations  on  this  process,  for  example  in  which  the 
resin  sets  chemically  during  the  injection  process  rather  than  during  a  separate  thermal  stage.  As 
a  new  manufacturing  process,  there  are  many  unresolved  technical  issues,  a  number  of  which  are 
under  active  study.  Representative  issues  include  the  strength  and  durability  of  the  finished  part, 
residual  stresses  and  warpage  resulting  from  thermal  gradients  induced  during  cure,  and  the  design 
of  the  fill  process  itself.  Within  the  latter,  important  process  variables  include  the  fill  rate  {i.e.  the 
driving  pressure  gradient),  the  capillary  number,  the  resin-fiber  wettability,  the  contact  angle,  the 
fiber  mat  characteristics,  the  uniformity  of  the  fiber  fill,  especially  at  corners  and  edges,  and  the 
location  of  the  inlet/outlet  ports  and  their  driving  pressure  histories. 

Among  the  important  issues  these  variables  seek  to  control  is  the  possible  formation  of  large  areas 
not  contacted  by  resin  (dry  spots).  Dry  spots  have  a  catastrophic  influence  on  strength  and  result 
in  an  unacceptable  product.  Also  very  significant  is  the  control  of  small  voids,  or  microbubbles, 
that  can  reduce  strength  by  a  factor  of  two. 

The  mathematical  analysis  of  RTM  is  described  at  the  continuum  level  by  Darcy’s  law 

V  =  — AVP  ,  Vv  —  0  , 

where  v  is  the  velocity  of  the  resin  as  it  is  injected  through  the  fiber  mat,  P  is  the  driving  pressure, 
and  X  =  K/fi  IS  the  transmissibility,  with  K  permeability  and  fa  resin  viscosity.  The  injection 
is  usually  modeled  as  a  free  boundary  problem,  with  the  resin  front  moving  at  the  above  Darcy 
velocity. 

A  number  of  authors  have  studied  the  microflow  properties  of  resin  moving  around  and  between 
fibers  twisted  into  a  string,  or  “tow”,  around  the  tows  woven  into  a  mat,  and  between  the  mats 
pressed  in  layers  into  the  mold  form  [Astrom  et  al.,  1992,  Bayramli  and  Powell,  1992,  Chen,  1993]. 
In  general  there  has  been  little  study  connecting  such  microflow  issues  to  continuum  level  descrip¬ 
tions  such  as  a  Darcy’s  law  flow  description  or  a  predictive  manufacturing  design  capability  to  select 
or  control  the  various  process  variables  for  optimization  of  manufacturing  design.  This  is  precisely 
the  issue  we  have  addressed. 

The  investigators  and  coworkers,  in  collaboration  with  a  manufacturing  research  team  from 
Northrup-Grumman.  initiated  such  a  study,  with  the  specific  goal  of  identifying  the  process  variables 
that  would  allow  reduction  of  microvoids.  A  preliminary  study  was  completed  [Chui  et  al.,  1995a]. 
The  study  used  a  new  continuum  level  description  of  the  flow  process.  Following  the  ideas  of 
porous  media  flow,  a  relative  as  well  as  an  absolute  permeability  was  introduced,  to  model  the  flow 
properties  of  the  microbubbles.  Thus  we  consider  the  Buckley-Leverett  equation  [Lake,  1989] 


St  +  V  -  vf{s)  ^  0  . 
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Here  .s  is  the  volume  fraction  of  the  void  space  occupied  by  the  injected  resin,  v  is  the  velocity 
from  Darcy’s  law  above,  and  /  is  a  fractional  flow  function,  that  expresses  resin  flow  as  a  fraction 
of  total  (resin  and  air)  flow,  /  is  expressed  directly  in  terms  of  the  relative  permeabilities,  ki{s), 
i  —  a  (air)  or  z  =  r  (resin),  and 


f  = 


kr  /  fJ>r 


ka/ f^a  k-p / fip 


while  A  above  is  replaced  by 


A  = 


We  use  a  standard  quadratic  form  kr  —  5“  for  the  resin  relative  permeability,  but  the  air  relative 
permeability,  ka  =  {sr  —  if  s  <  sr  and  ka  =  0  otherwise,  is  a  quadratic  as  a  function  of  the 
reduced  saturation  defined  in  terms  of  a  residual  air  saturation  sr.  We  take  sr  to  be  pressure 
dependent  to  express  a  higher  mobility  of  the  microbubbles  at  increased  pressure,  to  allow  for 
reduced  volume  of  the  (compressible)  bubbles  in  a  higher  pressure  field. 

The  resulting  flow  simulations  [Chui  et  ah,  1995a,  Chui  et  ah,  1995b]  were  shown  to  agree  quan¬ 
titatively  with  experiments  performed  in  Sweden  [Lundstrom  and  Gebart,  1992],  see  Figure  1,4. 


Void  Distribution  at  Breakthrough 
Outlet  Pressure:  1  atm 


Void  Fraction  (%) 


Normalized  Distance  to  Outlet 


Figure  1.4.  In  the  experiments  the  mold  was  filled  till  “breakthrough”,  after  which 
time  the  resin  cured  and  the  void  distribution  was  measured  along  the  mold.  The 
voids  concentrate  on  the  outlet  side  of  the  mold.  Our  relative  permeability  models 
have  a  few  physical  parameters.  The  figure  shows  that  these  parameters  can  be 
adjusted  to  give  quantitative  agreement  between  experiment  and  simulation. 


ADV.  COMP.  Sz  MODELING  FOR  NONLINEAR  STOCHASTIC  WAVES 


11 


The  most  important  flow  variable  for  the  control  of  microvoids  in  these  experiments  was  the 
evacuation  level  of  the  mold  during  filling.  Similar  experiments  were  performed  at  Northrup- 
Grumman,  reaching  the  same  conclusion.  It  was  found  that  components  of  excellent  quality  were 
produced  at  sufficient  level  of  evacuation  of  the  mold  prior  to  the  resin  injection. 


1.4.3.  Elastic  Plastic  Deformations.  The  proposers  and  coworkers  have  suggested  a  new  approach 
for  the  computational  modeling  of  elastic  plastic  deformation  [Glimm  et  ah,  1993b].  This  effort  is 
intended  to  allow  improved  simulations  for  large  deformation  problems,  as  occur  in  penetration 
mechanics,  tool  cutting,  punching,  and  other  metal  forming  processes.  The  new  methods  proposed 
are  useful  either  individually  or  in  combination.  The  principle  components  of  our  approach  are  (a) 
a  fully  conservative,  Eulerian,  upwind  finite  difference  approach  with  (b)  material  models  bcised  on 
a  hyperelastic  strain  energy,  (c)  front  tracking,  (d)  improved  physics  modeling  for  rate  independent 
plasticity,  and  (e)  special  modeling  to  allow  tracking  of  shear  bands  [Glimm  et  al.,  1995b].  The 
fully  conservative  Eulerian  formulation  [Plohr  and  Sharp,  1989,  Plohr  and  Sharp,  1992]  has  been 
completed,  and  represents  the  first  occasion  in  which  fully  conservative  equations  have  been  writ¬ 
ten  for  plastic  deformation  (or  even  for  the  Eulerian  formulation  of  elastic  deformation).  In  the 
case  of  plasticity,  the  equations  contain  new  physical  principles,  as  the  conservative  formulation 
implies  jump  conditions  not  contained  in  nonconservative  formulations.  These  two  formulations 
are  equivalent  for  smooth  flows  only.  The  conservative  independent  variables  are  strain  quantities, 
as  opposed  to  the  nonconservative  stresses  usually  used  as  primitive  variables.  The  formulation  is 
Eulerian,  so  that  Lagrangian  mesh  tangling  will  not  arise,  even  in  the  case  of  large  deformations. 
In  some  Lagrangian  codes,  mesh  tangling  is  avoided  by  use  of  the  method  of  erosion.  In  the  erosion 
method,  mesh  cells  are  simply  removed  from  the  computation,  together  with  the  material  contained 
in  them,  when  they  become  too  tangled  or  too  distorted  to  allow  the  computation  to  continue.  The 
approach  we  propose  will  use  mesh  and  interface  constructions  consistent  with  recognized  physical 
principles. 

The  fully  conservative  formulation  requires  a  hyperelastic  equation  of  state,  in  which  material 
properties  such  as  stress  are  represented  as  partial  derivatives  of  a  thermodynamic  strain  energy. 
The  hyperelastic  formulation  is  not  conventional  in  the  materials  modeling  community,  in  contrast 
to  incremental,  or  hyperelastic  thermodynamics,  in  which  thermodynamic  quantities  are  obtained 
from  the  solution  of  differential  equations.  In  contrast,  the  modeling  of  fluids  has  been  based  al¬ 
most  exclusively  on  a  hyperelastic,  or  strain  energy  based  equation  of  state  for  several  decades.  To 
formulate  the  hyperelastic  strain  energy,  we  follow  conventional  ideas,  and  use  a  small  deviatoric 
shear  model  [Garaizar,  1989]  in  which  the  isotropic,  or  fluid  equation  of  state  is  obtained  from  con¬ 
ventional  fluid  ideas,  such  as  a  stiffened  gamma  gas  law  or  a  tabular  equation  of  state  [ANO,  ]  and 
the  deviatoric  shear  strain  energy  is  added  perturbatively,  as  a  quadratic  correction  to  the  isotropic 
energy.  In  this  we  follow  conventional  ideas,  cf  Steinberg  and  Lund  [Steinberg  and  Lund,  1989]. 
Validation  is  required,  due  to  the  change  from  a  hyperelastic  (nonconservative,  incremental)  equa¬ 
tion  of  state  model  to  a  hyperelastic  one  [Wang  et  ah,  1993a]. 

Some  of  the  material  constants  for  rate  dependent  plasticity  can  be  obtained  by  direct  measure¬ 
ment,  while  others  are  inferred  only  indirectly  by  the  agreement  between  one  dimensional  plate 
impact  computations  and  experiments.  In  the  process  of  the  validation  study,  we  found  a  defi¬ 
ciency  in  the  original  physics  model  and  its  numerical  implementation,  which  implicitly  set  the 
plasticity  to  be  rate  independent  for  strain  rates  exceeding  the  Peierls  threshold.  This  physics 
modeling  assumption,  which  is  not  consistent  with  general  ideas  of  physics,  was  not  stated  clearly 
as  a  condition  of  the  model,  but  was  an  incidental  aspect  its  implementation.  This  model  seems  to 
be  replicated  in  several  plasticity  production  codes.  An  improved  physics  model  was  proposed,  and 
agreement  between  the  conservative  vs.  nonconservative,  and  hyperelastic  vs.  hyperelastic  formu¬ 
lations  was  obtained  [Grove  et  al.,  1993b.  Wang  et  al.,  1995].  Agreement  between  simulation  codes 
based  on  two  different  plasticity  formulations  means  that  the  material  parameters  that  allowed 
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[cm/|Lis]  Test  161  yi  [^^1/^5]  (b)  Test  232  yi  [cm/|is]  Test  390 


Figure  1.5.  Comparison  of  the  experimental  data  and  our  calculations  with  the 
new  material  model  with  the  new  material  parameter  C2  =  0.035  Mb  ms.  The  plots 
display  the  velocity  of  the  right-hand  free  surface  vs.  time  in  tests  161,  232,  and 
390. 

the  simulations  to  reproduce  the  experiments  were  the  same  for  the  two  formulations.  For  our 
improved  model,  one  of  the  indirectly  inferred  material  constants  required  a  value  different  from 
that  used  previously.  Agreement  between  simulation  and  experiment  was  satisfied,  both  before  and 
after  improvement  of  the  physics  model.  See  Figure  1.5  for  the  agreement  between  simulation  and 
experiment,  for  the  corrected  theory  and  the  corrected  material  parameter. 

In  order  to  compute  plate  impact  problems  in  an  Eulerian  formulation,  material  boundaries 
not  aligned  with  mesh  boundaries  were  encountered.  To  handle  this  difficulty,  a  simple  form  of 
front  tracking  was  introduced  [Wang  et  al.,  1993b].  In  order  to  extend  tracking  to  other  waves  of 
interest  in  plastic  deformation  problems,  a  model  of  shear  bands,  with  jump  relations  across  the 
shear  band  is  required.  We  have  begun  a  modeling  program  to  study  shear  bands.  Shear  bands 
arise  as  an  instability  in  laminar  shear  deformation.  As  one  portion  of  the  shear  layer  suffers  more 
plastic  deformation,  it  becomes  hotter,  hence  softer.  Thus  further  shear  intensifies  in  this  region, 
and  an  instability  breaking  the  laminar  behavior  is  started.  Shear  bands  in  metals  are  very  thin 
1  micron),  and  hence  impossible  to  resolve  as  part  of  a  macroscopic  deformation  simulation. 
They  must  be  modeled.  Our  model,  as  presently  completed  [Glimm  et  al,  1993b],  concerns  a  fully 
formed  shear  band  in  its  inner  region,  in  which  heat  conduction  is  in  approximate  balance  with 
heat  formation.  A  simplified  model  equation  for  the  evolution  of  the  shear  band  is 

6t  =  k6yy  =  cr7p 

where  6  is  temperature,  k  is  thermal  conductivity,  a  is  the  dominant  component  of  the  shear  stress, 
7  is  the  shear  strain  rate,  and  the  subscript  p  denotes  the  plastic  component.  The  key  aspect  of 
this  equation  is  the  nonlinear  source  term  describing  the  plastic  flow  rule,  z.e.,  the  heat  generated 
by  plastic  work.  Our  expression 

7p  =  ^sgn(a)ia|V-5(0)-'/- 

for  this  term  contains  some  modeling  assumptions.  See  [Clifton  et  al.,  1984,  Wright  and  Batra,  1985, 
Wright  and  Walter,  1987,  Walter,  1992]  for  related  studies. 

1.4.4.  Object  Oriented  N'limericol  Methods.  The  increasing  power  and  size  of  the  front  tracking 
code  has  led  to  a  corresponding  increase  in  the  time  and  effort  needed  to  train  new  personnel 
in  the  use  and  development  of  the  code.  Joint  research  efforts  with  the  Los  Alamos  National 
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Laboratory  and  the  Army  Research  Laboratory  also  require  the  addition  of  new  capabilities  in 
the  code  such  as  molecular  mixing,  elasto-plasticity,  and  radiation  hydrodynamics.  These  and 
other  upgrades  require  an  improved  implementation  of  the  front  tracking  code  using  modern  object 
oriented  programming  techniques.  This  work  has  led  to  improved  modularity  in  the  code  structure 
that  allows  new  applications  to  be  inserted  more  easily.  Recent  work  includes  the  installation  of 
elasto-plasticity  models  and  multiple  fluid  mixing  capabilities  in  the  front  tracking  code. 

2.  Professional  Activities 


2.1.  James  Glimm. 

2.1.1.  Editorial  boards. 

•  SIAM  Review  1989-1997 

•  Acta  Mathematicae  Applicandae 

•  Progress  in  Mathematical  Physics 

•  Applied  Mathematics  Letters 

•  Proceedings  American  Mathematics  Society 

•  Zeitschrift  fur  Angewandte  Mathematik  und  Physik 

•  Journal  of  Nonlinear  Analysis 

•  Computational  Geosciences 

•  Electronic  Research  Announcements  of  the  American  Mathematical  Society 

•  What’s  Happening  in  the  Mathematical  Sciences 

•  Scientific  Advisor:  Los  Alamos  National  Laboratory  Preprint  Bulletin  Board,  Computational 
Methods/Data  Analysis/ Wavelets/Lattice  Gases 

•  Series  Editor:  Mathematical  Methods  in  the  Life  Sciences,  World  Scientific  Publishers,  Inc. 

2.1.2.  Committee  Memberships  (Selected). 

•  Board  on  Mathematical  Sciences,  National  Research  Council  (1988-1994) 

•  Science  Policy  Committee,  SIAM  (1985-1996);  Chairman  1988-92 

•  Program  Director,  SIAM  Activity  Group  on  Geosciences 

•  Board  of  Timstees,  MSRI  (Berkeley)  (1989-1994) 

•  Advisory  Board,  AHPCRC  (Minnesota)  (1991-) 

•  Public  Information  Resource  Committee,  Joint  Policy  Board  on  Mathematics  (1987-  ) 

2.1.3.  Contacts  with  Industry.  Professor  Glimm  has  initiated  and  maintained  numerous  active  con¬ 
tacts  with  industry.  The  include  a  Resin  Transfer  Model  project  with  Northrup-Grumman,  reser¬ 
voir  modeling  with  the  Chevron  Corporation,  semiconductor  manufacture  with  the  DawnTech 
Corporation,  and  parallel  computing  projects  with  the  Intel  Corporation.  Industrial  and  dual  use 
technologies  research  began  with  a  survey  of  potential  uses  of  mathematical  modeling  in  advanced 
manufacturing  and  management  practices.  One  application  was  the  problem  of  geometric  pattern 
recognition  for  the  study  of  binding  sites  for  the  design  of  pharmaceuticals.  Here  the  large  combi¬ 
natorial  complexity  of  the  pattern  matching  optimization  function  is  the  main  difficulty.  A  rapid 
algorithm  was  developed  for  the  solution  of  this  problem.  The  method  was  compared  to  X-ray 
crystallography  to  determine  the  degree  to  which  hydrogen  bonds  alone  provide  binding  specificity 
for  binding  of  protein  on  DNA. 

2.1.4.  Contact  with  Army  Laboratory  Personnel.  Contacts  with  Drs.  Tim  Wright  and  John  Walters 
at  the  ARL  have  led  to  the  establishment  of  the  front  tracking  code  used  by  researchers  at  Stony 
Brook  at  the  laboratory.  Ongoing  collaboration  in  the  modeling  of  elasto-plastic  flows  and  adiabatic 
shear  bands  are  being  conducted. 
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Gerald  E.  Cooper  chief  of  the  Advanced  Research  Projects  Office  U.S.  Army  Concepts  Analysis 
Agency  has  initiated  contact  with  Professor  Glimm  on  the  modernization  and  optimization  of  the 
GEM  theater  combat  model.  Investigations  for  possible  areas  of  collaboration  are  ongoing. 

2.2.  John  W.  Grove. 


2.2.1.  Editorial  boards. 

•  Computers  and  Mathematics  with  Applications 

2.2.2.  Conferences  Organized. 

1.  Fifth  International  Conference  on  Hyperbolic  Problems,  Theory,  Numerics,  Applications,  Uni¬ 
versity  at  Stony  Brook,  Stony  Brook  NY,  June  12-17,  1994 

2.  Conference  in  Honor  of  the  Sixtieth  Birthday  of  James  Glimm,  University  at  Stony  Brook, 
Stony  Brook  NY,  April  20-22,  1995 

3.  Fifth  International  Workshop  on  the  Physics  of  Compressible  Turbulent  Mixing,  University 
at  Stony  Brook,  Stony  Brook  NY,  June  12-17,  1995 

2.2.3.  Collaboration  with  Government  Laboratories.  During  the  period  of  this  grant  Grove  made 
several  visits  to  government  laboratories,  both  to  the  Los  Alamos  National  Laboratory,  and  the 
Army  Research  laboratory.  Grove’s  work  at  these  sights  concerns  collaboration  with  laboratory 
researchers  (R.  Benjamin,  K.  Lackner,  R.  Menikoff,  and  D.  H.  Sharp  of  LANL,  J.  Walter  and 
T.  Wright  of  ARL)  on  numerical  methods  for  modeling  complex  flows.  In  particular  his  work  at 
Los  Alamos  is  primarily  based  on  using  front  tracking  to  model  complex  multi-fluid  mixing,  and 
my  work  at  ARL  is  devoted  to  the  development  of  front  tracking  methods  for  modeling  elasto- 
plasticity.  Both  collaborations  with  DOE  and  ARL  personnel  have  been  established  on  a  long  term 
and  continuing  basis.  The  details  of  these  research  projects  are  described  above. 

2.2.4.  Contacts  with  Industry.  The  INTEL  corporation  has  agreed  to  donate  eight  nodes  for  their 
Paragon  parallel  supercomputer  for  use  in  front  tracking  simulations  of  complex  mixing. 

2.2.5.  National  Service. 

1.  Member  National  Science  Foundation  Review  Panel  for  Small  Business  Innovation  Research 
Program,  September  26,  27  1994 

2.  Member  National  Science  Foundation  Review  Panel  Early  Career  Development  Program, 
February  21-23.  1995 

2.2.6.  Educational  and  Human  Resources  Development.  During  the  academic  year  of  1994/95  Grove 
supervised  four  graduate  students,  B.  Boston,  M.  J.  Graham,  A.  Pass,  and  S.  Tael,  two  of  whom 
are  women.  B.  Boston  will  complete  his  Ph.D.  at  the  end  of  June  1995,  while  A.  Pass,  and  S.  Tael 
expect  to  complete  their  Ph.D.’s  by  the  end  of  August  1995.  B.  Boston  has  accepted  a  position 
as  a  postdoctoral  associate  at  the  University  of  Stony  Brook  for  the  academic  year  of  1996-97,  A. 
Pass  has  accepted  a  position  at  Equitable  Insurance  Company  as  an  actuary,  S.  Tael  has  applied 
for  a  postdoctoral  position  at  the  Army  Research  Laboratory.  During  the  same  period  Grove  co¬ 
supervised  five  postdoctoral  students,  R.  Holmes,  X.  Lin,  K.-M.  Shyue,  R.  Young,  and  Y.  Zeng. 
R.  Holmes,  who  received  his  Ph.D.  in  August  1994  under  Grove’s  supervision,  has  been  awarded 
an  NSF  postdoctoral  fellowship  and  will  work  at  New  York  University  during  the  academic  year 
of  1995-96.  K.-M.  Shyue  received  a  full  time  position  at  the  University  of  Taiwan,  and  X.  Lin,  R. 
Young,  and  Y.  Zeng  will  continue  their  postdoctoral  positions  at  Stony  Brook  for  the  next  academic 
year. 
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3.  Publications  and  Technical  Reports 

ll]  B.  Boston.  J.  W.  Grove,  L.  F.  Henderson.  R.  Holmes,  D.  H.  Sharp,  Y.  Yang,  and  Q.  Zhang. 
Shock  induced  surface  instabilities  and  nonlinear  wave  interactions.  Report  No.  SUNYSB- 
AMS-93-20.  State  Univ.  of  New  York  at  Stony  Brook,  1993.  Proceedings  of  Eleventh  Army 
Conference  on  Applied  Mathematics  and  Computing. 

[2]  Y.  Deng,  J.  Glimm,  and  D.  H.  Sharp.  Mixing  and  chaotic  microstructure.  Los  Alamos  Science 
21:124-132,  1993. 

13]  J.  Grove,  R.  Holmes,  D.  H.  Sharp.  Y.  Yang,  and  Q.  Zhang.  Quantitative  theory  of  Richtmyer- 
Meshkov  instability.  Phys.  Rev.  Lett..  71(21):3473-3476,  November  1993. 

[4]  J.  Grove,  B.  Plohr,  D.  Sharp,  and  F.  Wang.  A  Godunov  scheme  for  elasto-plasticity.  In  Trans¬ 
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pages  113-130.  Army  Research  Office.  1993. 
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